#============================================================================
# 终极修复版：累积存活/死亡曲线（匹配箱线图/拟合曲线风格）
# 核心修复：1.解决dplyr select参数错误 2.A/B字母标注 3.全闭合边框 4.大轴数字 5.短刻度线
# 依赖：ggplot2 dplyr readr tidyr gridExtra grid
#============================================================================
library(ggplot2)
library(dplyr)
library(readr)
library(tidyr)
library(gridExtra)
library(grid)

# -------------------------- 1. 基础配置（匹配箱线图/拟合曲线风格） --------------------------
# 1.1 读取数据
data <- read_csv("data.csv", show_col_types = FALSE)

# 1.2 统一配色（与生命表/箱线图一致）
year_cols <- c("#1b9e77", "#d95f02", "#7570b3", "#E41A1C")
names(year_cols) <- c("total", "2015", "2020", "2025")

# 1.3 核心主题（关键调整：全闭合边框+大轴数字+短刻度线）
theme_matched <- theme_bw() +
  theme(
    # 基础字体（匹配箱线图，轴数字字号增大）
    text = element_text(family = "serif"),
    axis.title = element_text(size = 14, face = "bold", family = "serif"),  # 轴标题字号不变
    axis.text = element_text(size = 13, family = "serif", color = "black"), # 轴数字字号增大（关键调整）
    plot.title = element_text(size = 12, face = "bold", family = "serif"),
    # 关键调整1：全闭合边框（四周完整黑色边框）
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
    plot.background = element_rect(color = NA, fill = "white"),
    # 关键调整2：缩短刻度线（长度0.15cm，匹配箱线图/拟合曲线）
    axis.ticks = element_line(color = "black", linewidth = 0.6),
    axis.ticks.length = unit(0.15, "cm"),
    # 图例（匹配箱线图风格）
    legend.title = element_text(size = 12, face = "bold", family = "serif"),
    legend.text = element_text(size = 11, family = "serif"),
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.margin = margin(2,2,2,2),
    # 其他匹配项
    panel.grid = element_blank(),
    plot.margin = unit(c(0.8, 0.8, 0.8, 0.8), "cm")
  )

# -------------------------- 2. 生命表构建（与主程序一致） --------------------------
build_life_table <- function(df, class_col, name) {
  # 改用基础R过滤，避免dplyr语法兼容问题
  df_clean <- df[!is.na(df[[class_col]]) & is.numeric(df[[class_col]]) & df[[class_col]] > 0, ]
  counts_table <- table(df_clean[[class_col]])
  counts <- data.frame(
    age_group = as.numeric(names(counts_table)),
    nx = as.integer(counts_table)
  )
  counts <- counts[order(counts$age_group), ]  # 基础R排序
  
  # 核心指标计算
  total_n <- sum(counts$nx)
  counts$lx <- counts$nx / total_n * 1000
  counts$dx <- c(counts$lx[-nrow(counts)] - counts$lx[-1], counts$lx[nrow(counts)])
  counts$qx <- ifelse(counts$lx == 0, NA, counts$dx / counts$lx)
  counts$kx <- c(log(counts$lx[-nrow(counts)]) - log(counts$lx[-1]), 0)
  counts$method <- name
  
  return(counts[, c("age_group", "nx", "lx", "dx", "qx", "kx", "method")])
}

# -------------------------- 3. 生成生命表数据（BDT/IVC） --------------------------
# BDT（株丛径）
lt_diam_total <- build_life_table(data, "diam_class_year", "total")
lt_diam_year <- do.call(rbind, lapply(2015:2025, function(y) {
  if(y %in% unique(data$year)) build_life_table(data[data$year==y, ], "diam_class_year", as.character(y))
}))
lt_diam <- rbind(lt_diam_total, lt_diam_year)
# 基础R过滤，避免dplyr filter兼容问题
lt_diam <- lt_diam[!is.na(lt_diam$lx) & !is.na(lt_diam$qx) & !is.na(lt_diam$kx), ]
lt_diam$method <- factor(lt_diam$method, levels = c("total","2015","2020","2025"))
lt_diam$group <- "BDT"

# IVC（重要值）
lt_imp_total <- build_life_table(data, "imp_class_year", "total")
lt_imp_year <- do.call(rbind, lapply(2015:2025, function(y) {
  if(y %in% unique(data$year)) build_life_table(data[data$year==y, ], "imp_class_year", as.character(y))
}))
lt_imp <- rbind(lt_imp_total, lt_imp_year)
# 基础R过滤，避免dplyr filter兼容问题
lt_imp <- lt_imp[!is.na(lt_imp$lx) & !is.na(lt_imp$qx) & !is.na(lt_imp$kx), ]
lt_imp$method <- factor(lt_imp$method, levels = c("total","2015","2020","2025"))
lt_imp$group <- "IVC"

# -------------------------- 4. 计算累积存活/死亡指标（改用基础R） --------------------------
calculate_cumulative <- function(life_data) {
  # 按method分组计算，改用基础R实现
  methods <- unique(life_data$method)
  result_list <- list()
  
  for(m in methods) {
    sub_data <- life_data[life_data$method == m, ]
    sub_data <- sub_data[order(sub_data$age_group), ]  # 按龄级排序
    
    # 计算S_x（阶段存活率）
    sub_data$S_x <- NA
    if(nrow(sub_data) > 1) {
      for(i in 1:(nrow(sub_data)-1)) {
        if(!is.na(sub_data$lx[i+1]) & sub_data$lx[i] != 0) {
          sub_data$S_x[i] = sub_data$lx[i+1] / sub_data$lx[i]
        }
      }
    }
    
    # 计算S_i（累积存活率）
    sub_data$S_i <- 1  # 首项为1
    if(nrow(sub_data) > 1) {
      for(i in 2:nrow(sub_data)) {
        if(!is.na(sub_data$S_x[i-1])) {
          sub_data$S_i[i] = sub_data$S_i[i-1] * sub_data$S_x[i-1]
        } else {
          sub_data$S_i[i] = sub_data$S_i[i-1]
        }
      }
    }
    
    # 计算F_i（累积死亡率）
    sub_data$F_i <- 1 - sub_data$S_i
    # 数值修正（0-1之间）
    sub_data$S_i <- pmin(pmax(sub_data$S_i, 0), 1)
    sub_data$F_i <- pmin(pmax(sub_data$F_i, 0), 1)
    
    result_list[[m]] <- sub_data
  }
  
  # 合并结果
  result <- do.call(rbind, result_list)
  rownames(result) <- NULL  # 重置行名
  return(result)
}

# 计算累积指标
lt_diam_cum <- calculate_cumulative(lt_diam)
lt_imp_cum <- calculate_cumulative(lt_imp)

# 保存数据（基础R合并）
life_all_cum <- rbind(
  cbind(lt_diam_cum, type = "diam"),
  cbind(lt_imp_cum, type = "imp")
)
write.csv(life_all_cum, "FIG03_LifeTable_With_Cumulative.csv", row.names = FALSE)

# -------------------------- 5. 绘制累积曲线（最终定稿样式） --------------------------
# 5.1 通用绘图函数（关键修复：改用基础R数据处理，避免dplyr select错误）
plot_cum_final <- function(data, group_name, letter) {
  # 基础R选择列，避免dplyr select错误
  plot_data <- data[, c("age_group", "method", "S_i", "F_i")]
  # 移除NA值
  plot_data <- plot_data[!is.na(plot_data$S_i) & !is.na(plot_data$F_i), ]
  
  # 转换为长格式（基础R实现，避免tidyr兼容问题）
  plot_data_long <- data.frame(
    age_group = rep(plot_data$age_group, 2),
    method = rep(plot_data$method, 2),
    indicator = c(rep("Cumulative survival", nrow(plot_data)), rep("Cumulative mortality", nrow(plot_data))),
    value = c(plot_data$S_i, plot_data$F_i)
  )
  # 过滤NA值
  plot_data_long <- plot_data_long[!is.na(plot_data_long$value), ]
  
  # 绘图（保留所有样式要求）
  p <- ggplot(plot_data_long, aes(x=age_group, y=value, color=method, linetype=indicator)) +
    # 线条/散点（匹配箱线图/拟合曲线）
    geom_line(linewidth = 1, alpha = 0.8) +
    geom_point(size = 1, alpha = 0.8) +
    # 配色/线型
    scale_color_manual(values = year_cols, name = "Year") +
    scale_linetype_manual(values = c("solid", "dashed"), name = "Index") +
    # 坐标轴（匹配箱线图刻度密度）
    scale_x_continuous(breaks = seq(1, max(plot_data_long$age_group), 1), expand = c(0.02, 0)) +
    scale_y_continuous(expand = c(0.02, 0), limits = c(0, 1.05), breaks = seq(0, 1, 0.2)) +
    # 标签（A/B字母作为副标题，匹配箱线图）
    labs(
      x = "Age class", 
      y = "Probability", 
      title = group_name,
      subtitle = paste0("(", letter, ")")
    ) +
    # 主题（全闭合边框+大轴数字+短刻度线）
    theme_matched +
    # 调整副标题（A/B字母）位置：左上角，加粗大号
    theme(
      plot.subtitle = element_text(size = 16, face = "bold", family = "serif", 
                                   hjust = 0, vjust = 1, margin = margin(-20, 0, 0, 10)),
      plot.title = element_text(hjust = 0.5, margin = margin(0, 0, 10, 0))
    )
  
  return(p)
}

# 5.2 绘制BDT(A)/IVC(B)曲线（最终样式）
p_bdt_final <- plot_cum_final(lt_diam_cum, "BDT", "A")
p_ivc_final <- plot_cum_final(lt_imp_cum, "IVC", "B")

# 5.3 组合图表（匹配箱线图布局比例）
p_combined_final <- grid.arrange(
  p_bdt_final, p_ivc_final,
  ncol = 2,
  nrow = 1,
  widths = c(1, 1),
  heights = c(8),
  # 总标题（匹配箱线图样式）
  top = textGrob("Cumulative Survival and Mortality", 
                 gp = gpar(fontsize = 16, fontface = "bold", family = "serif"),
                 y = 1, vjust = 0)
)

# 5.4 保存图表（匹配箱线图/拟合曲线输出参数）
ggsave("FIG03_Cumulative_Curves_Final.png", p_combined_final, 
       width = 14, height = 8, dpi = 600, bg = "white")
ggsave("FIG03_Cumulative_Curves_Final.pdf", p_combined_final, 
       width = 14, height = 8, device = cairo_pdf)

# -------------------------- 6. 输出结果 --------------------------
cat("\n=== 累积曲线终极修复版绘制完成 ===\n")
cat("★ 核心修复项：\n")
cat("  1. 解决dplyr select参数错误（改用基础R数据处理）\n")
cat("  2. A/B字母标注：左上角括号格式（(A)/(B)），16号加粗serif字体\n")
cat("  3. 全闭合边框：图表四周完整黑色边框，线宽0.8\n")
cat("  4. 坐标轴数字：字号增大至13号（轴标题保持14号）\n")
cat("  5. 刻度线：缩短至0.15cm（匹配箱线图/拟合曲线长度）\n")
cat("★ 输出文件：\n")
cat("  - 累积指标生命表：FIG03_LifeTable_With_Cumulative.csv\n")
cat("  - 最终版累积曲线图：FIG03_Cumulative_Curves_Final.png/pdf\n")